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Abstract. The thermal emission from isolated neutron stars is not well understood. The X-ray spectrum is very close to a 
blackbody but there is a systematic optical excess flux with respect to the extrapolation to low energy of the best blackbody 
fit. This fact, in combination with the observed pulsations in the X-ray flux, can be explained by anisotropies in the surface 
temperature distribution. We study the thermal emission from neutron stars with strong magnetic fields B > 10'^ G in order to 
explain the origin of the anisotropy. We find (numerically) stationary solutions in axial symmetry of the heat transport equations 
in the neutron star crust and the condensed envelope. The anisotropy in the conductivity tensor is included consistently. The 
presence of magnetic fields of the expected strength leads to anisotropy in the surface temperature. Models with toroidal 
components similar to or larger than the poloidal field reproduce qualitatively the observed spectral properties and variability 
of isolated neutron stars. Our models also predict spectral features at energies between 0.2 and 0.6 keV for B = 10'^ - 10'"*. 
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1. Introduction 

Neutron stars (NS) with large magnetic fields (B > lO'^ G), 
the so-called magnetars, are becoming more and more abun- 
dant as new observations reveal phenomena that can only 
be explained by the action of strong magnetic fields. It is 
now believed that the small population (4 objects) of soft 
gamma repeaters (SGRs) are young neutron stars with mag- 
netic fields in the range x 10'"* - 10'^ G. Another subclass 
of candidates to be magnetars are the anomalous X-ray pul- 
sars (AXPs), whose high X-ray luminosities and fast spindown 
rates make them different from isolated radio pulsars or from 
NS in accreting X-ray binaries. The six rn embers of this family 
jTiengo et al. 200^ iMcGarrv et al. 20051) exhibit spin periods 
in the range 5-12 s, and their inferred magnetic fields (from 
their period derivative) are in the same range as SGRs (see e.g. 
IWoods & Thompson 20051 for a comprehensive review about 
these two families of magnetar candidates). 

A third rare family of NS, the radio-quiet isolated neutron 
stars among which RX J 1856.4 -3754 is the first and bright- 
est example ( Walter et al. 1996), shares some common fea- 
tures with the standard magnetars (SGR, AXPs): periods clus- 
tered in the range 5-10 s., and increasing evidence of large 
magnetic fields (> 10'^ G). The properties of the seven con- 
firmed members of this family are summarized in Tabled The 
most puzzling feature is the apparent optical excess flux (com- 
pared to the extrapolation of the best fit to the X-ray emis- 
sion) observed in several objects, which needs of the exis- 
tence of large temperature variations o yer the surface t o rec- 
oncile the optical and X-ray spectra /Pons et aL 2002h . The 
evidence of anisotropic temperature distribution is also sup- 



ported by the fact that several of the thermal spectra show clear 
X-ray pulsations with pulsation amplitudes from 5 to 20 %, 
while others (RX J1856, RX J1605) have upper limits of 1.3- 
3% to the maximum puls ation amplitude jBurwitz et al. 200*31 
van Kerkwiik et al. ^004*), which can be explained in terms of 
different relative orientations between the rotation and mag- 
netic field axis. Thus, it is worth to investigate the influ- 
ence of strong magnetic field configurations on the tempera- 
ture distribution, which has been shown to be able to creat e 
large anisotropies in neutron star crusts ( iGeppert et al. 20041) . 
or in the envelope, where further complications due to 
quantizing effects of the magneti c field or accreted mate- 
rial have been studied in detail jPotekhin et al. 20031) . But 
there is yet another important issue regarding the thermal 
emission on magnetized neutron stars. Below some critical 
temperature (depending on the composition and the mag- 
netic field strength), the gaseous layers of highly magne- 
tized neutron stars may undergo a pha se transit ion that 
turns the gas into liquid or solid state llLai 200 ll) . which 
strongly reduces the emissivity from the NS surf ace compared 
to the blackbody case (Brinkmann 1980: Turolla et al. 20041: 
IPerez-Azorfn et al. 2005[ivan Adelsberg et al. 20051) . 

In this paper our aim is to extend previous works on the 
anisotropies and thermal emission of magnetized neutron stars 
( Geppert et al, 200 4) in two main ways: by extending to lower 
density the calculations, within the model of a condensed sur- 
face, and exploring the effect of toroidal components of the 
magnetic field. The generation of toroidal fields in the early 
stages of a NS life, and its interplay with the poloidal com- 
ponent is a complex problem linked to convect ive instabili- 
ties, turbulent mean-field dynamo tBonanno et al. 2003.) or the 
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Hall instability jRheinhardt et al. 20041 ). The magnitude of the 
toroidal fields is unknown but usually thought to be larger than 
the poloidal component and, as we discuss in this paper, have 
interesting observational implications. 

This paper is organized as follows: in the next section the 
plasma properties in magnetic neutron stars are reviewed. In 
section 3, the magnetic field configurations used in the calcu- 
lations are described. In section 4, we describe the equations 
governing the thermal evolution and structure in the presence 
of large magnetic fields, the numerical code used for the calcu- 
lations and some tests. The microphysics input is discussed in 
section 5 and finally, in section 6, we present our results. 

2. Plasma properties in magnetic neutron stars 

The transport properties and in general all physical properties 
of dense matter are strongly aff'ected by the presence of intense 
magnetic fields. Before discussing in details the microscopical 
properties of magnetic matter, we begin by reminding some 
typical definitions and quantities that serve as indicators of the 
relative importance of the magnetic field. 

The pressure in the crust and envelope is dominated by the 
contribution of the degenerate electrons. Consider an electron 
gas whose number density is tie- In the absence of magnetic 
field, the Fermi momentum pp, or equivalently the wave num- 
ber kp - pp/fi is 



1/3 



37r>Z 
Aniu 



1/3 



(1) 



where ot„ is the atomic mass unit, and we have assumed that 
the ions, with atomic number Z and atomic weight A, are com- 
pletely ionized. This assumption allows to relate He and the 
density p, p - '-jAniu. Defining the dimensionless quantity: 



xf^ — = 0.010066, ^ , 

nieC \ A I 



1/3 



(2) 



the Fermi energy is ef - m^c^ yjl + xj^ and the Fermi tem- 
perature is Tf - (ef - meC^)/kB — rrieC^i yjl + - 1). If the 
matter is at temperature T, the electrons are degenerate when 
T Tf. This condition is fulfilled in the whole NS except for 
the outermost parts. 

The magnetic field aff'ects the properties of all plasma com- 
ponents, specially the electron component. Motion of free elec- 
trons perpendicular to the magnetic field is quantized in Landau 
levels, which produces that the thermal and electrical conduc- 
tivities (as well as other quantities) exhibit quantum oscilla- 
tions. These oscillations change the properties of the degen- 
erate electron gas in the limit of strongly quantizing field in 
which almost all electrons populate the lowest Landau level. 
The electron cyclotron frequency corresponding to a magnetic 
field B is given by 



eB 



COB - 



(3) 



and the magnetic field will be considered strongly quantizing 
if the temperature of the electrons is T <sc and the density 



p < Pb, where 

1 E — — ~ 
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Here kB is the Boltzmann constant and — 
(eB/tic)^^^/iK^ V2) is the electron number density at which the 
Fermi energy reaches the lowest Landau level. The magnetic 
field is called weakly quantizing if T < Tb but p > pB- In this 
case the quantum oscillations are not very pronounced and 
occur around their classical value. The oscillations disappear 
for T » Tb 01 p » Pb and the field can be treated as classical. 

Let us turn now to the properties of the ions. In the absence 
of magnetic field, the physical state of the ions depends on the 
Coulomb parameter 



F: 



(Ze)2 0.23Z2/p\i/3 



kBTai 



T6 



(6) 



where a, - (3/47rn,y^^' is the ion-sphere radius, and pe and 
T(, are, respectively, the density and temperature in units of 10^ 
g/cm^ and 10'' K. When F < 1 the ions form a Boltzmann gas, 
when 1 < F < 175 their state is a coupled Coulomb liquid, 
and when F > 175 the liquid freezes into a Coulomb lattice. In 
general, the quantization of the ionic motion will be significant 
for temperatures lower than the Debye temperature, which is 
approximately (for ions arranged in a bcc lattice) 



0.45- 



3.5 X 10- 



(!)' 



and ojp^ is the ion plasma frequency 



47TZ^e^ni 



1/2 



(7) 



(8) 



In presence of strong magnetic fields, the electrons in an 
atom are confined to the lowest Landau level, the atoms are 
elongated and with larger binding energy and covalent bonding 
between them. Therefore, below some critical temperature (de- 
pending on the composition and the magnetic field strength), 
the gaseous layers of highly magnetized neutron stars may 
undergo a phase transition that turns the gas into liquid or 
solid state depending on the value of the Coulomb parameter 
F tiLai 200 1). For typical magnetic field strengths of 10'^ G, a 
Fe atmosphere will condensate for T < 0.1 keV while a H at- 
mosphere needs temperatures lower than 0.03 keV to undergo 
the phase transition to a condensed state. In such a condensed 
neutron star surface made of nuclei with atomic number Z and 
atomic weight A, the pressure vanishes at a finite density 



p, * 560 AZ-^'^B%^ gem 



(9) 



where B12 is the magnetic field in units of lO'^ G. 

In this latter case, matter is in solid state and phonons 
become an important agent to transport energy. When T > 
Tf,, many thermal phonons are excited in the lattice, and the 
phonons behave as a classical gas. However, if temperature is 
low, T < To, phonons behave as a Bose quantum gas and the 
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Table 1. Properties of isolated neutron stars observed by ROSAT, Chandra, and XMM-Newton dPons et al. 20021: iHaberl 20oi 
[Haberl et al. 2005; Kaplan et al. 2003; Kaplan & van Kerkwiik 2005). 



Source 


kT 


P 


P 


T 


Optical 


Optical 


Pulsation 


Eline 


Brfi/Bcjc 




(eV) 


(s) 


10-'- (s s-') 


10*yr 




excess factor 


amplitude 


(keV) 


10" G 


RXJ0420.0-5022 


45 


3.453 


<9 




B= 26.6 


< 12 


0.12 


0.329 


< 18/6.6 


RXJ0720.4-3125 


85 


8.391 


0.07 


0.6-2 


B= 26.6 


6 


0.11 


0.270 


2.4/5.2 


RXJ0806.4-4123 


96 


11.371 


< 2 




B> 24 




0.06 




< 14/? 


1RXSJ130848.6+212708/RBS1223 


95 


10.313 


< 6 




28.6 


< 5 


0.18 


0.3 


?/2-6 


RXJ1605.3+3249 


95 








B= 27.2 


11-14 


<0.03 


0.46 


?/9.5 


RXJ1856.4-3754 


60 






0.5 


V= 25.7 


5-7 


<0.02 






IRXS J214303.7+065419/RBS1774 


101 


9.437 






R> 23 




0.04 


0.70 


?/14 



number of thermal phonons is strongly reduced. Therefore, the 
Debye temperature To allows to discriminate the quantum be- 
haviour from the classical one. Another important parameter 
related to the phonon processes of scattering is the so-called 
Umklapp temperature, Tu = ToZ^^^e^ /(3tiVf), with being 
the Fermi velocity of electrons. 

All the previous properties and definitions are visualized 
and quantified in Fig.^ where we show, in a phase diagram for 
neutron star matter, the Fermi temperature (Tf) and Tb (at B = 
lO'^ G), the Debye temperature (To) and the Umklapp tem- 
perature. The central shaded band indicates the region where 
matter is in the liquid state, according to Eq. (|6jl and the con- 
dition 1 < r < 175. The region below that is the solid state 
(F > 175) and the region above that corresponds to the gaseous 
state (F < 1). For reference, we have included two realistic 
temperature profiles (dot-dashed lines), corresponding to two 
different core temperatures (10^ K and 10*^ K). The dashed 
lines appearing on the left-lower corner indicate the transitions 
to different ionization states (25, 20 and 15 free electrons per 
atom). The outer layers are composed of pure iron (Z - 26). 
For B = lO'"' G, Eq. Q gives ps = 5 x 10^ g/cm^, therefore, 
the field is strongly quantizing only at low densities, weakly 
quantizing in most of the envelope and classical in the crust 
(Tb < T). The zero pressure density as defined by Eq. is 
7 X lO'* g/cm^ 

3. Magnetic field structure. 

Although there is robust observational evidence that the exter- 
nal magnetic field is well represented by a dipolar configura- 
tion, the internal structure of the magnetic field in neutron stars 
is unknown, so that one has the freedom to prescribe arbitrary 
configurations. For weak magnetic fields (magnetic force neg- 
ligible relative to the pressure gradient or gravity) the deforma- 
tion of the star is very small and the particular field structure 
is not important. Finding consistent numerical solutions of the 
Einstein-Maxwell equations describing the structure of neutron 
stars endowed with a strong magnetic field, including the ef- 
fects of the Lorentz force and the curvature of the spacetime 
induced by the stress-energy tensor of the magnetic field, is a 
difficult problem only solved for purely poloidal configurations 
jBocauet et al. 199^ or ver y recently including t oroidal mag- 
netic fields as perturbations (loka & Sasaki 200?). In previous 
works it has been shown that to obtain a significant deformation 




102 



106 1010 
P (gr/cm') 



1012 1014 



Fig. 1. Phase diagram of neutron star matter. The solid lines 
show the characteristic temperatures Tf, Tb, To, and Tu, where 
Tf, Tb have been calculated for a magnetic field strength ofB = 
lO'^G. The dashed lines show the domains of partial ionization 
(at B - 0). The dot-dashed lines are two realistic temperature 
profiles for two different models with core temperatures of 10** 
K and 10^ K. 

of the star magnetic fields of the order of 10'^ G are required. In 
this work, for simplicity, and partially justified by the fact that 
most of our models will be force-free (in a Newtonian sense) 
and less strong (lO'"* - lO'^ G), we will consider a spherical 
neutron star described by a spherically symmetric metric of the 
form 



ds^ = -e^'^^'^dt^ + e^^^'^dr^ + r^(de^ + sin^ 0d^'^) 



(10) 



where e^^''^ is the lapse function, e'^ = (1 - 2m/ry^^^ is the 
space curvature factor, and we have taken G - c - 1. The 
usual equation of hydrostatic equilibrium of a relativistic self- 
gravitating fluid is 



dP _ ^^_^p^d<i) 
dr dr 



(11) 



where P is the pressure and p is the mass-energy density. 

Previous works about the effects of internal magnetic fields 
restrict the calculations to the simplest models (homogeneous, 
core dipole, purely radial in a thin layer), mainly to simplify the 
problem. We address to the review by Tsuruta (1998), (section 
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5) for an overview of the effects of the magnetic field on the 
thermal structure and evolution. 

However, it is known that neutron stars are born with dif- 
ferential rotation and that convective instabilit ies play a signif- 
icant role during the early stages of evolu tion ( iKeil et al. 1996t 
iMiralles et al. 2000tlMiralles et al. 20021) . Both diff^erential ro- 
tation and convective motions should lead to non-trivial 
magnetic field struc tures with non-zero toroidal components 
llBonanno et al. 2003I) . that will evolve according to 



dB „ . 
— = -V X —R ■ 

dt Utt 



(12) 



where we have used the Newtonian equations for simplicity. 
The relativistic versions of the induction equation for non- 
rotating and rotatin g neutro n stars can be found in the literature 
dCieppert et al. 2000: RezzoUa & Ahmedov 20041) . Above, R is 
the resistivity tensor. 

In the classical (non quantizing) relaxation time approxi- 
mation, the conductivities are related between them through the 
magnetization p arameter {iORTn) where m is the non-magnetic 
relaxation time ( Urpin & Ya kovlev 1980). Then, the magnetic 
field evolution equation (ea.ll2> can be written as follows 

with cr|| being the electrical conductivities parallel to the mag- 
netic field. The first term at the right hand side of the above 
equation describes Ohmic dissipation and the last term is the 
Hall-drift, which is not dissipative but affects the current con- 
figurations. In general, for strong magnetic fields, ojetq » 1, 
the Hall-drift cannot be neglected. Notice that 



C U)bTq 



AnBcrn Aneue 



(14) 



which does not depend on the relaxation time, making 
evident the non-dissipative character of the Hall term. 
Even if the initial magnetic field is purely poloidal, 
it will develop a toroidal part during the evolution 

/Naito & Koiima 1994^ 'Hollerbach & Riidiger 20041 

[Rheinhardt et al. 20 04; ICumming et al. 20 04i) and it is 
necessary to consider how it afi^ects the thermal structure 
properties of neutron stars. In the remaining of this section we 
describe the magnetic field configurations used in this work. 



3.1. Dipolar magnetic fields 

In some previous works, the structure of the magnetic field has 
been assumed to be poloidal, in which case the field can be 
conveniently described i n terms of the Stokes stream function 
dCennert & Urt^in 1 994t iMiralles et al. 1 998t IPage et al. 2000t 
lOeppert et al. 2004 ). In spherical coordinates, and writing the 
(p component of the vector potential as - S(r)sm6/r, 
the magnetic field components can be written in terms of the 
Stoke's function S(r) as follows 

Br - ^ cos 6 



B, ^ - 



sin 0dS(r,t) 
r dr 



(15) 



If we consider the static solution obtained by extending the 
vacuum solution to the center of the star we have 



Bf — BqR 



cos 6 

~3~ 



Bf, = 



BoR^ sine 



(16) 



which corresponds to the choice S(r) - BoR/r. This solution 
diverges at r = 0. For a magnetic field confined to the crust, 
S (r) should vanish in the core due to proton superconductivity. 
In this work we prefer not to restrict ourselves to use poloidal 
configurations and we use a more general structure as described 
in the next subsection. 

3.2. Force-free magnetic fields 

A difl'erent, less restrictive, way to prescribe the interior mag- 
netic field is to consider a family of force-free fields. A 
force-free field is the simplest model for the equilibrium 
magnetic fie ld in the s olar corona, above an active region 
of sunspots iILow 19 93*: Wiegelm ann & Neukirch 20031) . This 
class of fields are normally used to model the pre-flare coronal 
configurations and have the advantage to allow large fields and 
currents to exist simultaneously without exerting any force on 
the material. This helps to simplify the problem because it al- 
lows to use the spherical solution for hydrostatic equilibrium 
without magnetic fields and because this configurations are not 
subject to the Hall drift. 

A force-free magnetic field obeys 



B ■ V;u = 



(17) 
(18) 



where the second equation is a result of the first one and the 
Maxwell's equation V ■ B - 0. 

Let us consider an axially symmetric magnetic configura- 
tion. In spherical coordinates, with (0, (f>) being the angular co- 
ordinates with respect to the magnetic field axis of symmetry, a 
general magnetic field can be written in terms of the </> compo- 
nents of the potential vector and the magnetic field as follows: 



B 



1 5(sinM^) 1 d(rA^) 



rsinO 



86 



dr 



B, 



(19) 



Therefore, Eq. reads: 

sin6i5(sin0B^) \d{rB^) I ( d(rBg) 
~r do 



r dr r \ dr 
sin6i5(sin6'A^) I dirA^) 
r do ' r dr 



dB,. 
'do 



,B, 



(20) 



For simplicity, we will consider solutions with yU ^constant, so 
that Eq. jl8l l is automatically satisfied, although more general 
solutions exist. The equality of the r, 6 components in Eq. j20> 
is obviously satisfied if we take B^ - fiA^. By analogy with the 
core dipole, we try the ansatz A^ - sin 6 A{r), which leads to 
the following equation for the ^ components of Eq. ( I20> 



cfiA{r) 2 dA{r) 



dr^ 



dr 



+ 1/1 



-,\A(r) 



. 



(21) 
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This is a form of the Riccati-Bessel equation for I - 1, which 
has solutions of the form 

A(r) = a ji(x) + bni(x) (22) 

where a,b are constants, x - fir and ii(x) and ni(x) are spherical 
Bessel functions of the first and second kinds. For Z = 1 we 
have, explicitly: 



«i(x) 



sm X cos X 



x^ 

cosx 



X 

sinx 



The spherical Bessel functions of the first kind are regular in 
the origin { ji{x) oc x'), while the functions of the second kind 
diverge as n/(x) oc x"*'^''. 

Hence, a general interior solution that matches (continuity 
of the normal component of the magnetic field) with the vac- 
uum dipolar solution at the surface is 



B = C 2- 



. cos f 



-Air), 



sin 6* 5(rA(r)) 
r dr 



, /isin 9A{r) 



(24) 



where C - and Bn is the value of the magnetic field at the 

pole. This magnetic field can be obtained from the following 
potential vector 



A ^Cijur cos 6 A(r), 0, sin Air)) 



Notice that this general solution includes all simple con- 
figurations as particular limits. The core dipole can be recov- 
ered by taking the limit fi ^ and considering only the mi(x) 
functions (a = in Ea.l22>. Alternatively, in the limit ju — > 
but considering only the family of regular solutions ji{x), we 
arrive to the homogeneous magnetic field. In both cases, the 

component vanishes. We can also find solutions that match 
continuously with the two components of the exterior dipole by 
setting, / = 1, fl = cos(fiRs), and b - sm{fiRs). The family of 
solutions is parametrized by the value of /i, which can be inter- 
preted as a wavenumber. If we want to build crustal magnetic 
fields that match with an external dipole we only need to adjust 
the wavenumber to have a vanishing radial component in the 
crust-core interface (r - Rint) and continuity of B,- and Bg at 
the surface (r - Rs). These values of yu are the solutions of 



tan(ij(Rint- Rs)) ^ nRh 



(26) 



In Fig. |2]we show three examples of crustal magnetic fields 
for the first three solutions of Eq. ( I26> with Ri-ai -9.2 km and 
Rs = 12.247 km. 

In principle there are more solutions of Eq. ( I20> than the 
linear one (B^ = /M^) that we have adopted. A more general 
discussion about force-free configurations can be fou nd in the 
Hterature JLow 1993U Wiegelmann & Neukirch 2003*). 

Finally, we also want to mention that this general force-free 
solution can be readily extended to higher order multipoles (i.e. 
quadrupole). This can be done by replacing in Eq. J24> cos 9 
and sinfl by the corresponding Legendre polynomial and its 
derivative, and using the spherical Bessel functions of the same 
index I. Another advantage of this force-free solution is that, 
for constant electrical resistivity (although in a realistic NS this 



15 



(23) -10 



-15 






2 4 



6 8 10 12 14 
x(km) 



6 8 10 12 I 
«(kfn) 



6 8 10 
x(knn) 



Fig. 2. Projections of the field lines on the (r, 6) plane for force- 
free magnetic field configurations confined in the region be- 
tween Ri„, - 9.2 km and Rs = 12.247 km. From left to right 
fl - 0.577, 1.569,2.591km '. All configurations match contin- 
uously to an external dipole. 



approximation is not appropriate), one can readily estimate the 
diffusion time. The evolution equation (I13> is simplified to 



dB 



(25) = __B 



which solution reads 



(27) 



(28) 



Here, Bq is the initial magnetic field and the decay time is 



T'diff = 73^5 that for typical conditions in young neutron stars 
is lO*" - 10^ years. 

4. Thermal diffusion in highly magnetized neutron 
stars 

4.1. Equations. 

For very slowly rotating NSs, and neglecting the magnetic 
force (which vanishes for force-free fields) the thermal evolu- 
tion of neutron stars can still be described by the energy balance 
equation 



dT 



of 



(29) 



where C,, is the specific heat (per unit volume), F is the en- 
ergy flux and the source term (e) includes all energy losses and 
sources (neutrino emission, frictional or accretion heating, etc). 
The evolution equation can also be written in integral form ap- 
plying Gauss' theorem 



edV 



In the diffusion limit, the energy flux is given by 



(30) 



(31) 
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where k is the thermal conductivity tensor Defining a new vari- 
able f - e^^''^T, the components of the flux can be written as 
foUows 



Fr - -{Kr, 



dyt + —det) 
r 

d,.T + —dgf) 

r 



(32) 



where /c,y are the components of the thermal conductivity ten- 
sor. The (p component of the flux is not considered because we 
assume axial symmetry. 

In the presence of strong magnetic fields, the thermal con- 
ductivities are different in the directions parallel and perpen- 
dicular to the magnetic field. In the classical relaxation time 
approximation, and considering that only electrons carry heat, 
the ratio between the parallel and perpendicular conductivi- 
ties is related to the magnetization parameter (cobtq) as follows 
(Einjin & Yakovlev 198^ 



1 + (ajBTof 



(33) 



The heat conductivity tensor in spherical coordinates and with 
the polar axis coinciding with the axis of symmetry of the mag- 
netic field can be written as follows 



/ + (OJBTof 



brbg 
brbs 



brbg 

bgb^ 



brb^ 
bgb^ 
bl 





bg 



b^, 


-br 



-bg 
br 


(34) 



where / is the identity matrix and br, bg, b^ are the components 
of the unit vector in the direction of the magnetic field. 
With the above expression for k, the flux reads 

ei'('->f = -K^ [vf + (wbTo)2 [bVf)b + lobTo [b x Vf)](35) 

The Hall contribution to the heat flux is given by the last term 
on the right hand side of Eq. ( I35t . If the magnetic field geom- 
etry has only poloidal components, and the temperature dis- 
tribution does not depend on the az imuthal angle, d), th e di- 
vergence of the Hall term vanishes taeppert et al. 2004*) and 
it does not affect the energy balance equation M9\ . However, 
for a magnetic field structure with a toroidal component, this 
term contributes to the heat ffux, even in axial symmetry. In the 
following, to simplify notation, we will omit the tilde over the 
temperature and we will use the symbol T for the red-shifted 
temperature. 

4.2. Boundary conditions. 

Boundary conditions can be imposed in either the temperature 
or the ffux in the boundaries of our numerical domain. Only a 
few years after birth the inner core of a NS becomes isothermal, 
therefore, in the core-crust interface (r = /?,„,) we will impose 
a fixed core temperature (7^). At the surface we impose 



F(B, T, 0b) = a(B, T, 9b)o-T'^ 



(36) 



where cr is the the Stephan-Boltzmann constant and a{B, T, 0b) 
is the integrated emissivity that depends on the partic- 
ular model, and 6b is the angle between the magnetic 





Fig. 3. Flux factor a as defined in Eq. (I36t for a dipolar mag- 
netic field with Bp - 5 X lO'^G, with (left panel) and without 
(right panel) taking into account the effect of the motion of the 
ions (considered as free particles). 

field and the direction normal to the surface element. At 
the temperatures of interest and for magnetic fields in- 
tense enough to produce the condensation of the gaseous 
layers, the emissivity at low energies is strongly re- 
duced compared to t he blackbody cas e and depends on 
the orientation ang le (Brinkmann 1980; Turolla et al. 20041 
IPerez-Azorm et al. 2005: .van Adelsberg et al. 2005.) . 

In Fig. |3] we show the ffux factor with (left) and without 
(right) taking into account the effect of the motion of ions for a 
dipolar magnetic fiel d of fin = 5x lO'^G. Based on our previous 
results ( Perez- Azor m et al. 20051) . we have obtained a polyno- 
mial fit of a(B, T, 6b) as a function of T(, (temperature in units 
of 10^ K) and cos(6b) for different magnetic field strengths (rel- 
ative error < 2%) with the following form: 



'=1 j=i 



(Ob) 



(37) 



The fl,j coefficients for B = 10'^ G and 5 x 10'^ G with and 
without taking into account the effect of the motion of ions are 
presented in tables (|4j0. The case of isotropic emission (black 
body) can be recovered by setting a - 1. 

4.3. Numerical test. 

The numerical algorithm consists of a standard finite difference 
scheme fully implicit in time. The temperature is cell-centered, 
while the ffuxes are calculated at each cell-edge. In order to test 
the code, we have studied the evolution of a thermal pulse in 
an infinite medium (neglecting all general relativistic effects), 
embedded in a homogeneous magnetic field oriented along the 
z-axis. If the conductivity is constant in the medium, an an- 
alytical solution for the temperature profile at a time t is the 
following: 



exp 



4-Kj_t 



sin 6 + 



cos^ 



1 + (WbTo)2 



(38) 



where Tq is a constant (the central temperature at the initial 
time, ?()), Kj_ is the transverse conductivity and (cobTo) is the 
magnetization parameter. To check the accuracy of the method, 
we have compared the numerical evolution of the pulse with 
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total conductivity includes the contributions of three carriers, 

K-Ke+ Krad + Kph (39) 

where Kg is the electron conductivity, Krad is the radiative (pho- 
ton) conductivity and Kph is the phonon conductivity. In non 
magnetic neutron stars, heat is transported mainly by electrons 
in the crust and in the inner envelope and by photons near 
the surface, while the phonon transport is negligible. However, 
in the presence of strong magnetic fields, this situation may 
change. While for the transport along the magnetic field the 
phonon contribution is still negligible, in the transverse direc- 
tion the electron transport is drastically suppressed, and the 
phonon contribution may become the most important one. Let 
us consider each of this contributions separately. 



Fig. 4. Temperature profiles at different times comparing the 
analytic solution (solid) and the numerical evolution (stars) of 
a thermal pulse in a medium embedded in a homogeneous mag- 
netic field. For clarity, in the numerical solution we have shown 
only one out of every six grid points. The parameters of the 
simulation are (ObTq - 3 and k^ - 10^ (a.u), to = lO"'*. The left 
panel shows the evolution of polar profiles and the right panel 
corresponds to equatorial profiles. 

the analytical solution, for different values of the parameters 
Kj^ and cobTq. As boundary conditions, we prescribe the tem- 
perature corresponding to the analytical solution in the surface 
and we impose F = at the center. In Fig.l^we show the com- 
parison between the analytical (solid) and numerical (stars) so- 
lution and for a model with cjbTq = 3 and a:^ = 10^ (a.u). 
The grid resolution is 100x40 (radialxangular). The deviations 
from the analytical solution in all cases studied are less than 
0.1%. 

5. Microphysics input. 

The microphysical ingredients that enter in the transport equa- 
tions ( I29> and i31\ are the specific heat and the thermal con- 
ductivity. Strictly speaking, the specific heat is not needed to 
obtain stationary configurations, but we have chosen to evolve 
Eq. (I29> without sources with a fixed inner temperature until 
the stationary solution is reached. Therefore, by using realistic 
microphysics input we will obtain also information about the 
thermal relaxation timescales. 

The dominant contribution to the specific heat is that from 
electrons and ions. For electrons we use the formulae corre- 
sponding to a relativistic degenerate Fermi gas while for ions 
we follow van Riper (1991). The most important ingredient is, 
however, the thermal conductivity, which has contribution from 
electron, photon and phonon transport. In this section we sum- 
marize the expressions used in the simulations. 



5.1. Thermal conductivity. 

The region of interest covers a large range of densities, from 
the core-crust boundary 10^* g/cnP) to the surface, which is 
given by Eq. in the models of condensed atmosphere. The 



5.1 .1 . Electron transport 

In the crust and the envelope of a neutron star the trans- 
port properties are mainly determined by the process of 
electron scattering off strongly correlated ions. The study 
of the transport properties of Coulomb plasmas with and 
without magnetic fields has been a focus of attention for 
decades (e.g. Flowers & Itoh 1976; Urpin & Yakovlev 19^ot 
iKaminker & Yakovlev 1 98 it lltoh et al. 1984I) . For the enve- 
lope, we will use the expressions obtained by Potekhin (1999), 
who calculated the thermal and electrical conductivities of de- 
generate electrons in magnetized envelopes by means of an ef- 
fective scattering potential that takes into account multiphonon 
processes in Coulomb crystals and an appropriate structure fac- 
tor of ions in Coulomb liquids. For the crust, the practical ex- 
pression derived by Potekhin (1999) have been later general- 
ized ( Gned in et al. 2001) by including how the size and shape 
of nuclear charge affects the transport properties as well as re- 
considering the electron-phonon scattering processes. 

In our calculations we are using the results from Potekhin 
(1999), whose code is of public domain '. For pedagogical pur- 
poses, in order to make evident the effect of a large magnetiza- 
tion parameter, we now summarize the classical relations for 
degenerated electrons. Schematically, the thermal conductivity 
can be written in terms of some effective relaxation time, t^-, 
as follows 



Kij = 



3ef 



-Tijiep), 



(40) 



where T,y are interpreted as inverse effective colUsion frequen- 
cies. In the non-quantizing case, we can write explicitly the 
different components in terms of the magnetization parameter 

(WfiTo) 



To 



1 + (jJbTQ 



1 + ojbto 



(41) 



The three main electron scattering processes that play a role 
in our scenario are scattering off ions, electron-phonon scat- 
tering and scattering off impurities. Semi-analytic expressions 
and fitting formulae for the relaxation times and thermal con- 
ductivities along the magnetic field for all three processes were 

' vjww. ioffe . rssi .ru/astro/conduct/condmag .html 
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Fig. 5. Magnetization parameter (wbtq) against density for dif- 
ferent temperatures (from top to bottom 10'', 10^, 10^, 10** K) and 
B - lO'^G. The solid lines are calculated with an impurity 
parameter of Q - 0.1, the dashed lines are for homogeneous 
matter (Q - 0) and the dot-dashed lines correspond to highly 
inhomogeneous matter Q - \Q. 



derived by Potekhin & Yakovlev (1996). The total contribution 
of electrons to the thermal conductivity is then calculated as 



(42) 



In Fig.|5]we show the magnetization parameter, related to 
the anisotropy of the thermal conductivity by Eq. j33> . as a 
function of density and for different temperatures. For com- 
parison, we show results with different impurity concentration 
parameter Q = nimp(-Zimp - Z)^/ni. The dashed, solid, and dash- 
dotted lines correspond to 2 = 0, 0.1, and 10, respectively. For 
highly inhomogeneous matter (Q = 10), the magnetization pa- 
rameter is strongly reduced in the crust ((Dbtq ^ 1). 

At large temperature the total electron conductivity is 
weakly dependent on temperature. If the temperature drops be- 
low the Umklapp temperature (T <k Tu), the Umklapp pro- 
cesses are disallowed and k oc T^*. Therefore, at high tempera- 
ture the dominant process is the electron-phonon scattering but 
at low temperature the scattering off impurities becomes the 
dominant contribution. 

5.1 .2. Photon transport 

Radiative conduction becomes the most effective transport 
mechanism in the outermost layers of the envelope, where elec- 
trons are non degenerate. We employ the expressions derived 
by Potekhin & Yako vlev (2001) for fully ionized iron, who fit- 
ted previous results jSilant'ev & Yakovlev 198ril) . 

Free-free transitions and Thompson scattering off free elec- 
trons are the two contributions to the total radiative conductiv- 
ity, that is calculated according to 
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Fig. 6. Thermal conductivity due to electron transport (dashed 
line), phonon transport (dot dashed line), photon transport 
(double dot dashed) and total (continuous line). The left pan- 
els show the conductivity in the direction longitudinal to the 
magnetic field while the right panels show the transverse con- 
ductivity, which is strongly suppressed. 



For temperatures below 10'' K, the dominant contribution 
comes from free-free transitions, which scales as a; p^^T^^^. 
Notice that for T < 10^ K we have k± a; 2k\i. 

5.1.3. Phonon transport 

Energy transport by phonons is usually orders of magnitude 
less effective than the usual electron or radiative transport. 
However, in the situation that we are studying, the large 
anisotropy induced by the magnetic field can suppress electron 
thermal conduction in the perpendicular direction by factors of 
10^ - 10*. Under this circumstances, transport by phonons be- 
come important, since this processes will become the most ef- 
fective way to transport energy in the perpendicular direction. 
For this reason, we need to include it in our calculations. 

In a first approximation, we consider a very simplified 
model, in which the phonon distribution is characterized by a 
Debye spectrum and all the relaxation times are functions of 
the wave vector of one mode only. In this approximation, the 
lattice thermal conductivity can be expressed as 



Kph 



exp (x) 



(exp(x) - 1)^ 



-dx 



(44) 



(43) 



where c, is the sound speed. To is the Debye temperature, x 
is a dimensionless variable (x = hcij/kBT) and t is the com- 
bined relaxation time, whose reciprocal is the sum of the recip- 
rocal relaxation times for all scattering processes considered, 
Umklapp and impurity scattering processes (both dissipative) 



Perez-Azorin, Miralles & Pons: Anisotropic tliermal emission from magnetized NS 



9 



and the three phonon normal scattering which are non dissipa- 
tive (Holland 1963; Konstantinov et al. 2003.): 



^-1 _i_ -,--1 _i_ ^-1 

T„ + T, + T., 



(45) 



At temperatures T > To, the lattice conductivity is mainly 
determined by the Umklapp processes, and the integral (I44> can 
be approximated to the expression 



pc^a 

'at' 



(46) 



vI/3 



where a « {^f^) is the lattice constant. At lower tempera- 
tures (T < Td), dissipative processes make the conductivity to 
increase very rapidly and the inclusion of impurity and normal 
phonon scattering becomes necessary. These processes (which 
conserve the total momentum) cannot by themselves lead to 
a finite thermal conductivity, but do not allow very large heat 
currents to be carried by modes of long wavelength. In this low 
temperat ure limit the ther mal conductivity can be expressed in 
the form llCallawav 196ll) 



Kph 



i^Io"^^ exp (x)(exp (x) - 1) ^c/xj 

In^C, \ n j ^TolT ^^1^4 (;t;)(exp (x) - lY^d^ 



(47) 



where r^' 



Ziman limit (IZiman 19711) 



Kph 



we recover the 



(48) 



where Uo is the Debye frequency and Timp is a parameter that 
takes into account the diflFerent atomic masses of the impurities 
and the lattice deformation. 

According to this expressions, in the crust and inner enve- 
lope, the heat transport along the magnetic field is dominated 
by the electron component, while phonons become the main 
transport agent in the transverse direction (see Fig. |6}. If the 
core temperature is low enough (T < lO^K), the phonon con- 
tribution becomes very important also in the parallel direction, 
making the crust to be quasi isothermal. Near the surface, the 
radiative conductivity dominates in both directions. 

6. Results. 

Our aim is to find stationary solutions of the temperature dis- 
tribution in a given background magnetic field configuration. 
We assume that the inner core is isothermal, and that the diffu- 
sion time of the magnetic field (rdiff) is much longer than the 
relaxation time to reach thermal equilibrium so that the mag- 
netic field is kept fixed. We also assume that the sources or 
sinks of energy (if any) are effective only at longer timescales. 
This assumptions are justified because both the magnetic dif- 
fusion time (either Ohmic or ambipolar) and the cooling time 
are > 10^ years, while the typical time to achieve the stationary 
solution (starting from a constant temperature profile) is ^ 10^ 
years. Notice that the diffusion timescale when the Hall insta- 
bility occurs is about 10'* years, so that in this case one would 
need to consider the coupled evolution of the temperature and 



the magnetic fields. Instead of solving the equation V-F -Q di- 
rectly, we evolve Eq. ( I29> . without sources, until the stationary 
solution is reached. 

The main eff'ect of the magnetic field on the temperature 
distribution can be guessed by looking at the expression of 
the heat flux ( I35t . When the magnetization parameter is large 
(wfiTo » 1), the dominant contribution to the flux is propor- 
tional to {a)BTQ)^(B ■ Vr). Therefore, in order to reach the 
stationary configuration the temperature distribution must be 
such that the surfaces of constant temperature are practically 
aligned with the magnetic field lines {B ■ <s; 1). This is 
shown explicitly in the left panel of Fig. where we show 
the stationary solution for a purely poloidal configuration con- 
fined to the crust and the outer layers ( poloidal confined, PC 
in the following). This alignment is enforced in most of the 
crust and envelope, and only near the surface strong radial gra- 
dients are generated. When we introduce a toroidal component 
the situation changes, because the Hall term in Eq. ( I35t induces 
large meridional fluxes (order wbTq) which result in an almost 
isothermal crust. This is clearly seen in the central panel, that 
shows the temperature distribution for a force-free magnetic 
field (FF) with a toroidal component present in the outer layers 
(crust and envelope). For comparison, we also considered an- 
other non force-free model (right panel) which has a toroidal 
component confined to a thin crustal region (toroidal confined, 
TC in the following), with a maximum value of 2 x lO'^G. It 
acts as an insulator keeping a different temperature at both sides 
of the toroidal field. In the region external to the toroidal field, 
since only the poloidal component is present we see again the 
alignment of isothermal surfaces with the magnetic field lines, 
which would not happen if the component extend all the 
way up to the surface, as in the central panel. We must stress 
again that the poloidal field is the same in all three models 
(solid lines. Bp = lO'^ G), but the field lines have been omit- 
ted in the central panel for clarity. The dashed Unes are con- 
tours of constant B^. In the lower panel of Fig. Q we show the 
same results but stretching artificially the low density regions 
to make visible the gradients near the surface. A slight north- 
south asymmetry provoked by the Hall term is visible in the 
right panel. 

The core temperature for all models is 5x10^ K. Thus, the 
anisotropy induced by the field becomes important not only 
in the crust but also in the condensed envelope. The direct 
consequence is a non-uniform surface temperature distribution 
shown in Fig. |8j where we show the angular distribution of 
the surface temperature for several magnetic field configura- 
tions, all of them with the same surface magnetic field (dipolar. 
Bp = 10'^ G). For comparison, we have also included (thick 
solid line) the semi-analytic temperature distribution derived 
by Greenstein and Hartke (1983), 

= Jcos^ «B + ^ sin^ j (49) 

where 9b is the angle between the normal vector to the surface 
and the magnetic field. The figure compares the following mod- 
els: core dipole (dashed line), PC (dotted line), FF (thin solid 
line) and TC (dash-dot). Qualitatively, the purely poloidal con- 
figurations (core dipole, PC) look similar to the Greenstein & 




Fig. 7. Upper panel: Temperature distribution in the crust of neutron stars with different toroidal components. The poloidal 
component is the same in all models {Bp = 10'^ G) and it is confined to the crust (solid lines). The left panel shows results 
for a purely poloidal field, the central panel a force-free configuration, and the right panel corresponds to a toroidal component 
confined to a narrow region of the crust. In the two latter cases the dashed lines show contours of constant B^. The scale has been 
stretched about a factor 2 to enlarge the crustal region. Lower panel; Same as the upper panel but stretching the scale to enlarge 
the region of the envelope. 
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Fig. 8. Surface temperature profiles as a function of the polar 
angle for different magnetic field configurations with the same 
value at the pole Bp = 10'^ G. The core temperature for all 
models is 5x10^ K. The models considered are: core dipolar 
(dashed), PC (dotted), TC (dash-dot), and FF (thin soHd Hne). 
In all cases we have included phonon transport effects with 
Timp = 0.1. The temperature distribution of Greenstein and 
Hartke j49> is also shown for comparison (thick solid line). 
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Fig. 9. Surface temperature profiles as a function of the po- 
lar angle for the PC and TC configurations. The dashed lines 
have been calculated without taking into account the phonon 
contribution, while the solid lines correspond to models in- 
cluding phonon transport effects with an impurity parameter 
of Timp =0.1. 



Hartke solution of Eq. J49> . with quantitative differences of the 
order of 10-20%. The general structure (relatively large hot po- 
lar region, and narrow cool equatorial band) is reproduced by 
all models without toroidal components of the magnetic field. 
This situation changes when a toroidal magnetic field is in- 
cluded, as for example in the force free configuration. In mod- 
els with important toroidal components, the surface thermal 
distribution consists of a small hot polar region and a relatively 
large cooler (about a factor 2-3) area. The slight north-south 



asymmetry provoked by the Hall term is only visible in the TC 
model (compare thin solid line with dot-dashed line). 

In the models where large meridional gradients in the 
crustal region are not present the phonon contribution to the 
thermal conductivity is not relevant and varying the impurity 
concentration barely changes the results. This is not true for 
the most extreme models (PC, TC), where phonon transport 
can make a difference. In Fig.|5]we compare the surface tem- 
perature distribution in the PC and TC models when phonon 
transport is switched off. The solid lines coiTespond to models 
in which the phonon contribution to the thermal conductivity 
is included (we have taken Fj 



imp 



0.1), while the dashed lines 
show models obtained without including the phonon conduc- 
tivity. Despite the fact that the effect of phonons is evident in 
the TC model, we must point out that the general distribution 
(small hot polar cap, larger cooler region) remains similar. 

A simple way to understand the results presented in this 
section is based on the following arguments. The models can 
be generally classified in two subclasses: i) magnetic field con- 
figurations that result in almost isothermal crusts (core dipole, 
FF, homogeneous) and ii) configurations for which large crustal 
temperature gradients are present (PC, TC). The first subclass 
includes models without toroidal components but also models 
with toroidal components present in the whole crustal region 
(i.e. FF). As discussed at the beginning of this section, the Hall 
term in Eq. ( I35> is responsible of the meridional heat flux that 
smears out temperature anisotropics in the crust. For such mod- 
els, the surface temperature distribution is well reproduced by 
the classical Greenstein and Hartke formula ( I49> but noticing 
that the dependence of Qb with the polar angle Q is different for 
each model. For a core dipole, we have 



cos Qb - 



4cos-6» 



(50) 



(51) 



1 -t- 3 cos^ Q ' 
for a FF model 

2 4 cos^ Q 

cos Qn — — — — 

( 1 -H H- (3 - iP-R^) cos2 e 

and for a homogeneous magnetic field cos^ Qb - cos^ Q. We 
have checked that T - T(Qb) looks very similar in all three 
cases despite the apparent differences in the surface distribution 
T{Q). The size of the hot polar cap can be easily estimated for 
this models. If we define the angular size of the polar cap as the 
angle where the temperature has decreased a given factor (say 
a factor 2, for example) with respect to the polar temperature, 
we can solve for Q to obtain the polar cap size. 

Models that admit strong crustal temperature gradients (PC, 
TC) do not obey Eq. j49> , and in principle there is no simple 
way to obtain how the temperature varies with the polar angle. 
The only general rule is that a strong toroidal component is 
necessary to produce small hot polar caps. 

6.1. Effective temperature. 

In Fig.^Jwe show the dependence of the effective temperature 
on the core temperature and the magnetic field strength. The ef- 
fective temperature is defined as L = AnRlcrT^^, where L is the 
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Fig. 10. Dependence of the effective temperature on the core 
temperature and the magnetic field strength for two different 
configurations core dipole (solid) and force free (dashes). We 
show results for three different core temperatures, from bottom 
to top 10^ 5 X 10^ and 10'^ K. 

total integrated luminosity over the surface. This effective tem- 
perature is the quantity usually obtained from black-body fits 
to observational data, and plotted on cooling curves to compare 
data with theoretical predictions. The three solid lines corre- 
spond to three different core temperatures, from bottom to top 
10^, 5 X 10^, and 10** K, and for a core dipole configuration. 
The dashed lines correspond to the same core temperatures but 
for a force free magnetic field. In all cases we observe a sys- 
tematic lower effective temperature (a factor x 2) in configura- 
tions with toroidal magnetic fields in the crust-envelope region. 
This means that among NS with similar ages (i.e. similar core 
temperatures during the neutrino dominated cooling era), those 
with strong toroidal fields have an apparent effective temper- 
ature about a factor 2 smaller than those with low magnetic 
fields or purely dipolar configurations. 

Our result can also be compared to the classical formula 
that relates the temp erature at the base of the en velope with the 
surface temperature dOudmundsson et al. 1982h 

/ J, > 0.455 

n,8 = 1.288 (52) 

where T;, g is the temperature in the base of the envelope in 
units of 10*^ K, TsSfi is the surface temperature in units of 10'' 
K and gu is the gravity acceleration in units of 10''* c.g.s. For 
the three core temperatures (10^,5 x 10^, 10^ K) used in Fig. 
1101 and assuming an isothermal crust, the surface temperature 
of non-magnetized NSs would be of 0.05, 1.85 and 8.58 x 10^ 
K, respectively. The quite different effective temperatures pre- 
dicted by different models are relevant for the interpretation of 
the comparison of observational data with cooUng curves. 

6.2. Quantizing magnetic field effects. 

In the previous results the quantizing character of the mag- 
netic field has been neglected. For 10'^ G, Eq. (|5} gives ps ~ 



Fig. 11. Surface temperature profiles as a function of the polar 
angle for the same models as in Fig.|8lbut including quantizing 
effects on the electron thermal conductivity. 



4.8 X 10^ g/cm^, while the density of the condensed surface 
is ps » 7 X 10^ g/cm^*. Therefore only in the outermost thin 
layer (p < ps) the magnetic field can be considered strongly 
quantizing while in most of the envelope it is weakly quantiz- 
ing. Including the quantizing effects on the conductivities in 
a transport code is challenging because the rich structure in 
Landau levels makes necessary a robust code and high resolu- 
tion to handle properly the gradients that might develop near 
each transition. For a selected number of models we have per- 
formed the calculations including quantizing effects with the 
purpose of understanding the qualitative and quantitative dif- 
ferences with the classical case. In Figs. ^2 ™d [21 we show 
the results for the same models as in Figs. |8] and |9] but in- 
cluding quantizing effects on the electron thermal conductiv- 
ity (Potekhin, 1999). The two main facts that we observe in 
this figures (present as well in other models not shown) are the 
following. First, the average effective temperature is generally 
lower and the anisotropy is more pronounced, i.e. a smaller an- 
gular size of the hot polar region. Second, the surface temper- 
ature distribution shows small oscillations in models without 
toroidal components near the surface (core dipolar, PC, TC), 
associated to the oscillatory behaviour of the thermal conduc- 
tivity in the quantizing case. This can be explained by the fact 
that the poloidal component is practically radial and heat trans- 
port in the meridional direction is strongly suppressed. The 
radial gradients are different at each latitude due to the dif- 
ferent magnetic field strength, and therefore different densities 
at which electrons are filling the corresponding Landau levels. 
This is shown in Fig. [O] where radial temperature profiles for 
three different polar angles are plotted. The Landau levels are 
clearly visible. This oscillatory behaviour cannot be smoothed 
out by meridional heat fluxes because they are suppressed by 
a factor x (w^to)^. The exception is the FF model, for which 
oscillations are not observed, because the it has a toroidal com- 
ponent extended up to the surface. This makes possible the ex- 
istence of heat flux in the meridional direction because the Hall 
term is order {(jJbTo)k±. 
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Fig. 12. Same as Fig.|9lbut including quantizing effects on the 
electron thermal conductivity. 
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Fig. 13. Temperature profiles as a function of the density for 
different polar angles = 0, 45, 60, and 90°. The magnetic field 
configuration is TC with Bp - lO'^ G and a core temperature 
of 5x10^ K. 

6.3. Influence of the physical conditions of the outer 
layers. 

One of the important issues under debate is whether the enve- 
lope and atmosphere will be in a gaseous or condensed state. 
In order to estimate the dependence of our results on the choice 
of a particular model, we have compared four different outer 
boundary conditions which represent the different possibilities 
one can find. A common approach is to solve th e 2D heat 
transfer equation only in the crust (Geppert et al. 2004 ) and 
match at p ^ 10^" g/cm^ to sorn e magnetized envelope so- 
lution (" Potekhin & Yakovlev 200lV Alternatively, as we have 
discussed in this work, there is the condensed surface model 
in which the emissivity at low energies may vary depending 
on the way that the motion of ions i n the lattice (fixed or free 
ions are the two limits) i s treated dPerez-Azorm et al. 2005t 
Ivan Adelsberg et al. 20051) . Also one can consider the simplest 
model which is to assume that the condensed surface (e.g. 
Ps -7x10* g/cnv', for B = 10'^ G) radiates as a blackbody. We 



Fig. 14. Influence of the boundary condition on the surface 
temperature distribution for a model with Tc = 5 x lO^K 
and a force-free magnetic field configuration with B - 
lO'^G. Results are shown for blackbody emission (dot-dashed), 
gaseous magnetic envelope (dashed), and metallic surface with 
(solid) and without (triple dot - dash) taking into account the 
motion of the ions. 

have analyzed this four possibilities and we show the resulting 
surface temperature distribution for all four models in Figs. 1141 
and^l for the classical and quantizing cases, respectively. The 
main difference, as stated previously, is that including quan- 
tizing effects leads to lower average temperatures. Notice also 
that when quantizing effects are included, the size of the hot 
polar region is smaller and the temperature is nearly constant 
in a large part of the surface of the star. 

It must be stressed that these are all FF configurations, for 
which the crust is very close to isothermal, and the gradients of 
temperature are generated in the low density region. The con- 
clusion from this comparison is that not only the temperature at 
the base of the crust, but also the physical conditions in the low 
density layers affect the total luminosity (the average effective 
temperature). However, the general shape of the surface tem- 
perature distribution is qualitatively similar in all cases, which 
leads to conclude that irrespectively of the physical assump- 
tions, neutron stars with strong magnetic fields do have large 
surface temperature variations. 

It a forthcoming work we will study models with stronger 
magnetic fields (magnetar rather than isolated NS conditions) 
in which the quantizing effects are probably even more impor- 
tant. For the remaining of this work, and having established 
the qualitative trends that differentiate classical and quantiz- 
ing models, we will focus on analyzing, in the classical limit, 
a number of other different issues that might have important 
consequences on the emission properties. 

6.4. Influence of impurities. 

The influence that impurities or defects in the lattice have on 
the final temperature distribution may be important. Impurity 
scattering dominates either at very low temperatures (where 
phonon scattering is suppressed) or when the impurity level is 
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Fig. 15. Same as Fig.^^but including quantizing effects of the 
magnetic field. 



very high. For isolated NSs the values of the impurity parame- 
ter may vary from Q ^ 10"^ in very pure crusts io Q> 10 in the 
amorphous inner crust ( I Jones 20041) . In accreting neutron stars 
Q is set by the composition of the nuclear burning occur ring at 
low density, and it is Ukely that Q ^ 100 dSchatz et aTl 999). 
The impurity content also determines the critical field above 
whic h the Hall eff'ect dom inates over purely Ohmic dissipa- 
tion dCumming et al. 20041) . Given this uncertainty, we have 
explored a variety of models with the impurity parameter to 
test the sensitivity of our results to the impurity concentration. 

In Fig.[^we show the resulting surface temperature distri- 
butions corresponding to the PC and FF configurations and for 
different values of the impurity parameter Q - (dot-dashed), 
0.1 (dashed), and 10 (solid). Only small corrections to the PC 
configuration are visible, while the lines are indistinguishable 
for the FF model. Therefore, the exact value of the impurity 
concentration might be important for the long term evolution of 
the magnetic field and the temperature, but it does not seem to 
be crucial for the stationary solution corresponding to a back- 
ground field. 

6.5. Pulsations. 

The anisotropic temperature distribution obtained from our cal- 
culations will translate into periodic pulsations if the neutron 
star is rotating and the magnetic and rotation axis are not 
aligned. Within a fully relativistic framework that includes light 
bending effects (Page 1995; Page & Sarmiento 1996), we have 
calculated the visible luminosity curves for a number of mod- 
els with different magnetic field strengths and configurations. 
In Fig.[n]we show the observed luminosity obtained for mod- 
els with a core dipolar configuration with Bp = lO'"* G and 
Tc - 10^ K. We denote by O the angle between the observer 
and the rotation axis and by S the angle between the rotation 
and magnetic axis. The numbers next to each line are the max- 
imum pulsed fraction (MPF): 



MPF = 



F + F 

max ^ n: 



(53) 



Fig. 16. Influence of the impurity content for models with Tc = 
5 X lO^K and two different magnetic field configurations (PC 
and FF), both with the same poloidal component (B = 10'^ G). 
We show results for Q - (dash-dots), Q - O.l (dashes) and 
e = 10 (solid). 



The dependence of the MPF on the different parameters can 
be understood by analyzing the results summarized in Tables|2l 
(core dipolar configurations) and|3l(force free configurations). 
Several interesting conclusions can be drawn from the results. 
First, for a given core temperature, it depends very weakly on 
the strength of the magnetic field (Bp), but can by up to twice 
larger when the toroidal magnetic field is included (force free 
models). For poloidal fields the MPF in the models we ana- 
lyzed is 16%, but it can be increased to 25-30% in the force-free 
models. We must remind that the toroidal component is larger 
in about one order of magnitude than the poloidal one (see Eq. 
i24\ . fiRs ~ 10). The anisotropy and therefore the variability 
may be increased by using other magnetic field configurations 
with larger toroidal components. The observed variability of 
isolated neutron stars in consistent with this results, but some 
of them show large pulse fractions (1 1% in RX J0720, 12% in 
RX J0420, 18% in RBS 1223) than seem to indicate the exis- 
tence of toroidal interior magnetic fields, and large angles be- 
tween the rotation and magnetic axis. The lack of pulsations 
in (RX J 1856) is compatible with nearly aligned rotation and 
magnetic axis (S < 6°). The information about the variability 
correlated with the effective temperature and the optical excess 
flux can therefore give relevant information about the magnetic 
field structure. 



6.6. A comparison to blaMody models. 

Since most spectral fits to real data are made with simple 
blackbody models, we have taken one of our models (FF, 
Bp = 10'^ G, Tc = 5 X 10^ K) and fitted our results to 
a single blackbody. We have assumed a column density of 
Hfj - 1.5 X 10^° cm"^, typical of galactic interstellar medium 
absorption. The comparison between this model, and a BB 
fit is shown in Fig. The X-ray part of the spectrum is 
well fitted by a single blackbody but our model predicts an 
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Table 2. Maximum pulsed fraction for O - njl and S - n/l with a core dipolar configuration. 
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S„ = 3 X 10'^ G Bp = 10" G Bp = 2.5 x 10'^ G fi,, = 5 x 10'^^ G Bp = 10'* G 

Tc^lO^K 0.15 0.14 0.13 0.11 0.10 

re = 5xlO^K 0.16 0.16 0.16 0.15 0.13 

rc = 10'*K 0.16 0.16 0.16 0.16 0.15 



Table 3. Maximum pulsed fraction for O - njl and S - 7t/2 with a force-free configuration. 
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0.18 
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Fig. 17. Observed flux variability and the corresponding max- 
imum pulsed fraction for three different orientations of the 
magnetic axis with a fixed observer position {O - n/l): solid 
(S = 0), dotted (S = n/l) and dashed (S = 7r/4). The magnetic 
field structure is core dipolar with Bp = 10'^ G and Tc - 10' 
K. 



optical flux about a factor 4 larger than the blackbody fit to 
the high energy part. This factor may vary depending on the 
magnetic field strength and geometry and it is consistent with 
the systematic excess flux observed in the optical counter- 
parts of isolated neutron stars. More interestingly, the con- 
densed surface models also predict the existence of an edge 
at an energy E ^ h (cje^ + ayjJa>B^^ (van Adelsbera et al. 2005; 
IPerez-Azorm et al. 200.5h . where cjp = (Ane^ne/mey^^ is the 
electron plasma frequency, that for typical magnetic fields 
(10'3 - lOi'* G) faUs in the range 0.2-0.6 keV (depending also 
on the gravitational redshift). Some spectral features have been 
reported in that range although they are usually associated to 
proton synchrotron lines. The only object for which an inde- 
pendent estimate of the magnetic field is available is J0720, 
for which the measure of P - 6.98 ± 0.02 x lO'^'^js s' ^) im- 
plies B = 2.4 X 10'^ G (Kaplan & van Kerkwiik 20051) . The 
observed spectral feature is fitted by a Gaussian absorption 
line at an energy of 0.27 keV, and has been associated to cy- 
clotron resonance scattering of protons in a magnetic field with 
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Fig. 18. Comparison between the spectra of a FF model (solid 
fine) with Bp = lO'^G and T^ = 5 x 10' K and a single 
blackbody fit (dashed line). The parameters of the real NS are 
M - 1 AMq, R - 12.27 km, and we have taken d - 111 pc and 
riH - 1.5 X 10^"cm"^. We have assumed O - 7t/2 and S - 0. 
The parameters of the BB fit are given in the figure. The optical 
flux of the FF model is a factor 4.3 larger than the BB fit. 



B = 5 X 10'^ G (lHaherletal.2004h . Assuming a magnetic field 
strength (from P) of B - 2.4 x lO'^ G, the condensed sur- 
face model predicts a phase dependent edge at an energy (lo- 
cal) of 0.35 keV, which would imply a redshift of z = 0.29. The 
phase dependent emitted spectrum for one of our models (FF) 
is shown in Fig.[^ We have taken O - n/2 and S = 7r/2. The 
feature is strongly dependent on the orientation, being stronger 
when the magnetic field axis is pointing to the observer and 
practically undetectable when the magnetic axis is normal to 
the direction of observation. The angles have been chosen to 
show the most extreme case, where the variability is very large. 
As reported in Table |3] this particular model has a maximum 
pulsed fraction of 0.24. 

7. Conclusions. 

In this paper we have presented the results of detailed calcula- 
tions of the temperature distribution in the crust and condensed 
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Fig. 19. Phase dependent emitted spectrum (unabsorbed) of a 
FF model with Bp = lO'^G and = 5 x 10^ K. We have as- 
sumed O = n/2 and S = n/2. From bottom to top, the different 
lines correspond to phase angles of 0, 30, 45, 60, and 90°. The 
hot polar cap (and therefore the magnetic field axis) is pointing 
to the observer when the phase angle is 0° . 



envelopes of neutron stars in the presence of strong magnetic 
fields. The surface temperature distribution has been calculated 
by obtaining 2D stationary solutions of the heat diffusion equa- 
tion with anisotropic thermal conductivities. From the variety 
of strengths and configurations of the magnetic field explored, 
we conclude that variations in the surface temperature of fac- 
tors 2-10 are easily obtained with magnetic fields in the range 
{B > lO'^-lO'^ G). The average luminosity (or the inferred ef- 
fective temperature) does not depend much on the strength of 
the magnetic field, but it is drastically affected by the geome- 
try, in particular by the existence of a toroidal component. The 
toroidal field acts effectively as a heat insulator forcing heat to 
flow towards the poles. Therefore, it is the particular geome- 
try (a priori unknown) of the magnetic field what eventually 
determines the size of the hot polar caps. A back of the en- 
velope calculation to estimate the size of the polar cap is the 
following. For purely radial magnetic fields the non-magnetic 
solution is not affected while for purely tangential fields the 
temperature gradient is quite larger. Therefore the hot polar cap 
will be determined by the angular size of the region in which 
the magnetic field is nearly radial. For a classical dipole, the 
condition B^ = Bg leads to sin 6* « 2/ V5, which gives a hot po- 
lar cap of about 63°. For FF models, for example, the condition 
bI = bI + Bj leads to sin 6^2/ -jsV^jpW. The configurations 
we have employed correspond to fiR ^ 16, which gives an es- 
timate of the size of the hot polar cap of 7°, in good agreement 
with the numerical results. 

For purely poloidal configurations, the surface temperature 
is high in a large fraction of the star surface and lower in a 
narrow equatorial band, while for configurations with toroidal 
components of the same strength as the poloidal one the tem- 
perature distributions is more close to hot polar cap with a large 
cooler area at low latitudes. Thus, this latter family of models 
shows larger pulsation ampHtudes and optical excess flux, in 



very good agreement with the observed properties of isolated 
neutron stars. We defer to future work for detailed fitting of 
real data with our models, but preliminary calculations show 
that the spectral energy distribution and its variability can be 
easily explained without fine tunning of the model parameters. 
This can be interpreted as indirect evidence of the existence 
of toroidal fields in the crust and envelopes of NSs. We have 
also investigated the influence of some relevant inputs such as 
the physical conditions of the surface (condensed, gaseous) by 
varying the outer boundary conditions, i.e., the emissivity at a 
given temperature and B. We found that the main conclusions 
remain qualitatively unchanged, although quantitative differ- 
ences can arise. We have also explored the effect of having 
different impurity content, finding that their effect is not impor- 
tant in general, being only visible in models without toroidal 
components. 

Another interesting result is that the condensed surface 
models predict the existence of an edge at an energy E « 
H^ojbi +(jjp/ojB,^ that for typical magnetic fields falls in the 
range 0.2-0.6 keV, where some spectral features have been re- 
ported, and usually associated to proton synchrotron lines. The 
energy of the spectral feature observed in J0720, as well as 
its pulsation amplitude predicted by our models are consistent 
with the inferred magnetic field. We also plan to extend our 
work to calculations with stronger magnetic fields and higher 
temperatures, typical condition of magnetars (SGRs, AXPs). 
The mean caveat that we must point out is the large uncer- 
tainty in the particular structure of the magnetic fields inside 
neutron stars and the need of a full 2D calculation of the rela- 
tivistic structure of neutron stars with arbitrary magnetic fields. 
The bottom line is that magnetic fields do change significantly 
the thermal emission from isolated neutron stars and cannot be 
overlooked if one expects to infer valuable information (radius, 
gravitational redshift, composition) from the observed spectral 
energy distribution. What one infers from the blackbody fits to 
X-ray observations (assuming a known distance to the object) 
is the product T^Rl^, and model dependent variations in the es- 
timation of the effective temperature translate into the estimate 
of the radius. Our models with toroidal components result in in- 
ferred radii about a factor 3-5 larger than the BB radius or than 
the inferred radius from a model with only poloidal component. 
This naturally solves the problem of the apparent smaUness of 
some isolated neutron stars. 

If the existence of strong magnetic fields in isolated NSs 
is confirmed, we will need more detailed calculations coupling 
the evolution of the magnetic field with the temperature before 
we can establish firm constraints on NS properties by fitting 
observational data. 
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Table 4. Coefficients of the fit to the emissivity from a condensed surface (Ea.l37> for B - 5 x lO'^ G taking into account the 
effect of the motion of ions. 



Oij 1 2 3 4 5 6 



1 


0.451428 


0.652402 


-0.329039 


0.123172 


0.107814 


-0.0967540 


2 


-0.304007 


-0.626414 


-5.79256 


17.3578 


-20.8844 


9.57023 


3 


-0.378297 


4.48855 


3.42040 


-26.5822 


38.1295 


-18.8900 


4 


0.562032 


-3.96036 


-0.866190 


17.5741 


-26.3781 


13.3072 


5 


-0.226386 


1.36369 


0.0520957 


-5.29775 


8.03772 


-4.06952 


6 


0.0299993 


-0.169923 


0.0232908 


0.568170 


-0.880843 


0.450272 



Table 5. Same as table |4] without taking into account the effect of the motion of ions. 



Ojj 1 2 3 4 5 6 



1 


0.0520960 


0.906625 


-2.27950 


3.50014 


-2.73453 


0.831342 


2 


-0.092588 


3.96721 


-6.00109 


1.76058 


3.73908 


-2.65634 


3 


0.1884173 


-6.88457 


16.2169 


-20.2954 


13.1221 


-3.36603 


4 


-0.0905443 


5.4989 


-15.7136 


24.2137 


-19.2636 


6.14122 


5 


0.0197174 


-1.92340 


6.01936 


-9.94622 


8.30053 


-2.73754 


6 


-0.001500 


0.239160 


-0.782216 


1.33117 


-1.13048 


0.376769 



Table 6. Same as tableHwith B = 10'^ G. 





1 


2 


3 


4 


5 


6 


1 


0.37727 


-0.0297095 


0.172343 


0.934048 


-1.86330 


1.05495 


2 


-1.32853 


7.41110 


-16.7803 


15.8174 


-4.53211 


-1.57610 


3 


1.70549 


-7.45357 


16.4152 


-14.9966 


5.40906 


0.726946 


4 


-0.826587 


1.60020 


2.15163 


-13.0149 


14.5113 


-5.62832 


5 


0.185560 


0.229112 


-3.75523 


9.92179 


-9.52118 


3.30446 


6 


-0.0162908 


-0.076025 


0.652672 


-1.58083 


1.48421 


-0.505131 



Table 7. Same as table|5]with B = 10'^ G. 



fly 1 2 3 4 5 6 



1 


0.0682268 


0.945109 


-2.14347 


2.99701 


-2.20175 


0.647961 


2 


-0.269514 


4.23065 


-6.58908 


2.17098 


4.27706 


-3.28306 


3 


0.323785 


-3.69871 


2.05278 


6.96692 


-11.4661 


5.26209 


4 


0.030067 


-0.483140 


11.2452 


-27.8763 


26.7549 


-9.23072 


5 


-0.066218 


0.77907 


-6.42987 


14.4756 


-13.4195 


4.50075 


6 


0.0118741 


-0.131711 


0.949772 


-2.10127 


1.94095 


-0.648770 



